Point Judith 16S rRNA amplicons analysis
Rebecca Stevick
July 2023
About the Data
Conditions
3 oysters per bucket
3 tissue samples per oyster: gut, inner shell, outer shell
= 108 samples total
| Bucket | Location | Treatment | Station |
|---|---|---|---|
| 1 | Southern | Control | BHC |
| 2 | Southern | Control | BHC |
| 3 | Southern | Control | BHC |
| 4 | Southern | Enriched | BHC |
| 5 | Southern | Enriched | BHC |
| 6 | Southern | Enriched | BHC |
| 7 | Northern | Control | BC |
| 8 | Northern | Control | BC |
| 9 | Northern | Control | BC |
| 10 | Northern | Enriched | BC |
| 11 | Northern | Enriched | BC |
| 12 | Northern | Enriched | BC |
Sequencing
16S V6 region
Primers
1064R GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGCGACRRCCATGCANCACCT
967F-A TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCTAACCGANGAACCTYACC
967F-B TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCNACGCGAAGAACCTTANC
967F-C TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCAACGCGMARAACCTTACC
967F-D TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGATACGCGARGAACCTTACC
2 x 250 bp paired-end Illumina sequencing
Analysis
Initially done with QIIME2, dada2 denoising. Eukaryotic reads removed here.
Setup
Load libraries and settings
library(tidyverse)
library(scales)
library(ggpubr)
library(qiime2R)
library(vegan)
library(gt)
library(patchwork)
library(UpSetR)
library(gplots)
library(ggh4x)
library(ComplexHeatmap)
# set global theme
theme_set(theme_classic()+
theme(legend.position ="none",
panel.background = element_blank(), plot.background = element_blank(),
axis.title = element_text(size=18),
legend.title = element_text(size=18, face="bold"),
axis.text = element_text(size=15),
legend.background = element_blank(),
legend.text = element_text(size=18),
strip.text = element_text(size=22),
plot.title = element_text(hjust=0.5, size=22),
strip.background = element_rect(fill="grey80", color="transparent")))
set.seed(567)
# set global options for code output
knitr::opts_chunk$set(echo=TRUE, warning=FALSE,message=FALSE)Import data
# set file paths for metadata and qiime2 output
metadatafile <- "metadata/PJ_V6Samples_Metadata.txt"
tablefile <- "qiime2output/table.qza"
taxonomyfile <- "qiime2output/taxonomy.qza"
#treefile <- "qiime2output/rooted-tree.qza"
# make ASV table per sample
SVs<-read_qza(tablefile)$data
# make dataframe from metadata
metadata<-read_q2metadata(metadatafile) %>%
# clean up sampletype
mutate(SampleType=factor(SampleType, levels = c("gut","inner swab", "outer swab"),
labels = c("Gut", "Inner Shell", "Outer Shell")))
# make dataframe of taxomony per ASV
taxonomy<-read_qza(taxonomyfile)$data %>% parse_taxonomy()Clean up data
# vector of non-bacterial ASVs
eukASVs <- taxonomy %>%
filter(Kingdom=="d__Eukaryota" | Kingdom=="Unassigned" | Phylum=='Cyanobacteria' | Kingdom=='d__Archaea' ) %>%
rownames_to_column("ASVs")
# clean up ASV data
datafullASVs <-
# transpose ASV matrix and make dataframe
t(SVs) %>% as.data.frame() %>%
# set sample names as a column
rownames_to_column("SampleID") %>%
# add in the metadata
full_join(metadata) %>%
# make into long form
pivot_longer("4d42c811fce6eaabd10ae7ce7334bf7f":"1467c8a6122ac19b493fb587db2b32c9") %>%
# remove non-bacterial ASVs
filter(!name %in% eukASVs$ASVs)
# clean up ASV data
dataASVsClean <- datafullASVs %>%
# calculate percent abundance from counts
group_by(SampleID) %>% mutate(percent=value/sum(value)) %>%
# add in all taxonomy data
left_join(taxonomy %>% rownames_to_column("name")) %>%
# remove control data
filter(SampleType !="water" & SampleType!="control")
# aggregate data per level
taxasumsclean <- as.data.frame(SVs) %>%
# remove controls
select(-c(RS234_S134, RS235_S140, RS246_S110, RS247_S117)) %>%
rownames_to_column("ASV") %>%
filter(ASV %in% datafullASVs$name) %>%
column_to_rownames(var="ASV") %>% qiime2R::summarize_taxa(taxonomy)
# clean up ASV matrix table
SVsControls <- as.data.frame(SVs) %>%
# remove controls
select(-c(RS234_S134, RS235_S140, RS246_S110, RS247_S117)) %>%
# set ASV names as a column
rownames_to_column("ASV") %>%
filter(ASV %in% datafullASVs$name) %>%
# set ASV names as rownames, then transpose
column_to_rownames(var="ASV") %>% t()
# check that the ASVs in datafullASVscutoff match SVsControls
setdiff(colnames(SVsControls), unique(datafullASVs$name))## character(0)# Average per Bucket
SVsControlsAvg <-
as.data.frame(SVs) %>%
# remove controls
select(-c(RS234_S134, RS235_S140, RS246_S110, RS247_S117)) %>%
# set ASV names as a column
rownames_to_column("ASV") %>%
filter(ASV %in% datafullASVs$name) %>%
# convert reads to percent per sample
mutate(across(where(is.double), ~ .x/sum(.x))) %>%
# set ASV names as rownames, then transpose
column_to_rownames(var="ASV") %>% t() %>%
as.data.frame() %>% rownames_to_column("SampleID") %>% left_join(metadata) %>%
unite("BucketType", Bucket, SampleType) %>%
group_by(BucketType) %>%
summarise(across(`4d42c811fce6eaabd10ae7ce7334bf7f`:
`1467c8a6122ac19b493fb587db2b32c9`, mean)) %>%
column_to_rownames(var="BucketType") General QC
Sequencing controls
Mock expected
# expected mock abundances and plot
mockplot <- tibble(SampleName="Expected_Mock",value=0,
MockTaxa=factor(c("Pseudomonas aeruginosa", "Escherichia coli",
"Salmonella enterica", "Enterococcus faecalis",
"Lactobacillus fermentum","Staphylococcus aureus",
"Listeria monocytogenes", "Bacillus subtilis")),
percent=c(4.2,10.1,10.4,9.9,18.4,15.5,14.1,17.4)) %>%
ggplot(aes(x=SampleName, y=percent, fill=MockTaxa))+
geom_col(position="fill", alpha=0.8)+
theme(legend.text = element_text(size=12, colour="gray20", margin = margin(b = 10, unit = "pt")),
legend.position = "right",legend.direction = "vertical",
axis.ticks.y = element_line(inherit.blank=FALSE, color="grey30"))+
scale_fill_manual(values=c("#c5543a","#b8995a","#9cd250","#6eb386",
"#8c9fc3","#7e4abb","#c45891","#4b393e"))+
scale_y_continuous(labels = scales::percent_format(), expand=c(0,0))+
labs(y="Expected mock species abundance",x=NULL,fill="Mock species", title="Expected mock")
mockplot
Positive control
# sequenced blank and mock control plot
seqmockcontrol <-
# transpose ASV matrix and make dataframe
datafullASVs %>%
# select only positive control data
filter(SampleName=="MOCK.CON") %>%
filter(value!=0) %>%
# add in taxonomy and new column of ASV name
left_join(read_qza(taxonomyfile)$data, by=c("name"="Feature.ID")) %>%
unite("ASVname", c("name", "Taxon"), sep=": \n") %>%
# filter just top 20 taxa, group all Others
mutate(TaxaOther=forcats::fct_lump_n(f=ASVname, w=value, other_level="Others", n=10)) %>%
mutate(TaxaOther=reorder(TaxaOther, -value)) %>%
# calculate percentages per taxa per Sample
group_by(SampleName) %>% mutate(percent=value/sum(value)) %>%
# reorder based on expected mock
mutate(TaxaOther=factor(TaxaOther,
levels=c("281e3624d68f8ff9b997a85e9d5d9f56: \nd__Bacteria; p__Firmicutes; c__Bacilli; o__Bacillales; f__Bacillaceae; g__Bacillus",
"913c2544bf1019abcf31d51cac6e9edb: \nd__Bacteria; p__Proteobacteria; c__Gammaproteobacteria; o__Enterobacterales; f__Enterobacteriaceae",
"f0e13523941c32d798d606bd5df1d43e: \nd__Bacteria; p__Proteobacteria; c__Gammaproteobacteria; o__Enterobacterales; f__Enterobacteriaceae; g__Escherichia-Shigella",
"b8e412fddb124a840dc886233adb7289: \nd__Bacteria; p__Firmicutes; c__Bacilli; o__Lactobacillales; f__Lactobacillaceae; g__Lactobacillus",
"09a0e1e07fd7747e33cdfd76e592f3e6: \nd__Bacteria; p__Firmicutes; c__Bacilli; o__Lactobacillales",
"32ca4fa1d1056e4f185ab0220d55a757: \nd__Bacteria; p__Firmicutes; c__Bacilli; o__Lactobacillales; f__Listeriaceae; g__Listeria",
"144f9431e8a6187a98700a50e2a9af21: \nd__Bacteria; p__Proteobacteria; c__Gammaproteobacteria; o__Pseudomonadales; f__Pseudomonadaceae; g__Pseudomonas",
"446cc81456db29f56e72032cc90e8602: \nd__Bacteria; p__Proteobacteria; c__Gammaproteobacteria; o__Enterobacterales; f__Enterobacteriaceae",
"f9eed79f2bd69b0da2385144b13658a7: \nd__Bacteria; p__Firmicutes; c__Bacilli; o__Staphylococcales; f__Staphylococcaceae",
"5c9b89d5acf8e0020d9de19b53d21663: \nd__Bacteria; p__Firmicutes; c__Bacilli; o__Staphylococcales; f__Staphylococcaceae; g__Staphylococcus",
"Others")))
# make palette based on expected mock
pospalettess2<- seqmockcontrol %>%
mutate(colorpal = case_when(grepl("Escherichia",TaxaOther) ~"#9cd250",
grepl("446cc8",TaxaOther) ~ "#c45891",
grepl("Enterobacteriaceae", TaxaOther) ~ "#b8995a",
grepl("Pseudomonas",TaxaOther) ~ "#7e4abb",
grepl("Staphylococcaceae",TaxaOther) ~ "#4b393e",
grepl("Listeria",TaxaOther) ~"#8c9fc3",
grepl("Lactobacillales",TaxaOther) ~ "#6eb386",
grepl("Bacillus",TaxaOther) ~ "#c5543a",
TaxaOther=="Others" ~ "grey40"))
colors <- distinct(pospalettess2, TaxaOther, colorpal)
pal <- colors$colorpal
names(pal) <- colors$TaxaOther
posconperc <- seqmockcontrol %>%
ggplot(aes(x=SampleName, y=percent, fill=TaxaOther))+
geom_col(position="fill", alpha=0.8, color="white")+
theme(legend.text = element_text(size=12, colour="gray20", margin = margin(b = 10, unit = "pt")),
legend.position = "right",legend.direction = "vertical",
axis.ticks.y = element_line(inherit.blank=FALSE, color="grey30")) +
scale_fill_manual(values=pal) +
scale_y_continuous(labels = scales::percent_format(), expand=c(0,0))+
labs(y="Percent ASV abundance",x=NULL,fill="Positive control ASV", title="Positive")
posconperc
Negative control
negconperc <-
datafullASVs %>%
# select only negative control data
filter(SampleName=="NEG.CON") %>%
filter(value!=0) %>%
# add in taxonomy and new column of ASV name
left_join(read_qza(taxonomyfile)$data, by=c("name"="Feature.ID")) %>%
unite("ASVname", c("name", "Taxon"), sep=": \n") %>%
# filter just top 20 taxa, group all Others
mutate(TaxaOther=forcats::fct_lump_n(f=ASVname, w=value, other_level="Others", n=12)) %>%
mutate(TaxaOther=reorder(TaxaOther, -value)) %>%
# calculate percentages per taxa per Sample
group_by(SampleName) %>% mutate(percent=value/sum(value)) %>%
ggplot(aes(x=SampleName, y=percent, fill=TaxaOther))+
geom_col(position="fill", alpha=0.8)+
theme(legend.text = element_text(size=12, colour="gray20", margin = margin(b = 10, unit = "pt")),
legend.position = "right",legend.direction = "vertical",
axis.ticks.y = element_line(inherit.blank=FALSE, color="grey30"))+
scale_fill_manual(values=c("#A6CEE3", "#1F78B4", "#B2DF8A", "#33A02C", "#FB9A99", "#E31A1C", "#FDBF6F", "#FF7F00",
"#CAB2D6", "#6A3D9A", "#FFFF99", "#B15928","grey40"))+
scale_y_continuous(labels = scales::percent_format(), expand=c(0,0))+
labs(y="Percent ASV abundance",x=NULL, fill="Negative control ASV", title="Negative")
negconperc
Summary figure
mockplot+posconperc+negconperc+
plot_layout(guides="collect") & theme(legend.position = "bottom")
Reads per sample
reads <- dataASVsClean %>%
group_by(SampleID) %>% mutate(sumreads=sum(value)) %>%
# simplify data to one value per sample
distinct(Location, Treatment, SampleType, Group, SampleName, sumreads) %>%
# start plotting
ggplot(aes(x=SampleName, y=sumreads, fill=Group))+
facet_nested(.~SampleType+Location+Treatment, scales="free",space="free",
nest_line = element_line(color="white"),
strip = strip_nested(
text_x = elem_list_text(face = c("bold","bold","bold",NA,NA,NA,NA,NA,NA,
NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
size = c(24,24,24,20,20,20,20,20,20,
14,14,14,14,14,14,14,14,14,14,14,14))))+
geom_col(color="white")+
scale_fill_manual(values=c("lightgreen","darkgreen","lightblue","darkblue"),
labels=c("Northern Ambient", "Northern Enriched", "Southern Ambient", "Southern Enriched"))+
scale_y_continuous(expand=c(0,0), labels=scales::label_comma())+
labs(y="Number of sequences",x=NULL,fill=NULL)+
theme(axis.text.x = element_blank(), strip.text = element_text(size=14))
reads
Rarefaction curves
# number of ASVs per sample
(S <- specnumber(SVsControls))## RS126_S16 RS127_S28 RS128_S40 RS129_S52 RS130_S64 RS131_S76 RS132_S88
## 675 553 511 672 268 627 358
## RS133_S5 RS134_S17 RS135_S29 RS136_S41 RS137_S53 RS138_S65 RS139_S77
## 555 588 392 573 725 569 420
## RS140_S89 RS141_S6 RS142_S18 RS143_S30 RS144_S42 RS145_S54 RS146_S66
## 488 648 543 773 760 533 502
## RS147_S78 RS148_S90 RS149_S7 RS150_S19 RS151_S31 RS152_S43 RS153_S55
## 634 517 461 443 653 476 421
## RS154_S67 RS155_S79 RS156_S91 RS157_S8 RS158_S20 RS159_S32 RS160_S44
## 749 534 579 468 482 472 339
## RS161_S56 RS162_S68 RS163_S80 RS164_S92 RS165_S9 RS166_S21 RS167_S33
## 478 715 796 916 611 821 719
## RS168_S45 RS169_S57 RS170_S69 RS171_S81 RS172_S93 RS173_S10 RS174_S22
## 740 1048 770 702 722 861 608
## RS175_S34 RS176_S46 RS177_S58 RS178_S70 RS179_S82 RS180_S94 RS181_S11
## 789 863 867 964 811 796 643
## RS182_S23 RS183_S35 RS184_S47 RS185_S59 RS186_S71 RS187_S83 RS188_S95
## 608 731 785 786 712 503 722
## RS189_S12 RS190_S24 RS191_S36 RS192_S48 RS193_S60 RS194_S72 RS195_S84
## 829 688 704 850 795 640 569
## RS196_S96 RS197_S97 RS198_S104 RS199_S111 RS200_S118 RS201_S124 RS202_S130
## 685 880 1168 999 781 478 798
## RS203_S136 RS204_S142 RS205_S98 RS206_S105 RS207_S112 RS208_S119 RS209_S125
## 1465 1083 683 800 834 865 280
## RS210_S131 RS211_S137 RS212_S143 RS213_S99 RS214_S106 RS215_S113 RS216_S120
## 763 985 838 558 1054 573 709
## RS217_S126 RS218_S132 RS219_S138 RS220_S144 RS221_S100 RS222_S107 RS223_S114
## 784 846 582 904 894 780 757
## RS224_S121 RS225_S127 RS226_S133 RS227_S139 RS228_S145 RS229_S101 RS230_S108
## 771 283 1080 954 1045 915 940
## RS231_S115 RS232_S122 RS233_S128
## 796 453 795# smallest number of reads in a sample
(raremax <- min(rowSums(SVsControls)))## [1] 9490# rarefied number of taxa per sample
Srare <- vegan::rarefy(SVsControls, raremax)
# slope at the end of the rarefacetion curve per sample
Sslope <- vegan::rareslope(SVsControls, raremax)
# plot observed vs rarefied number of ASVs
plot(S, Srare, xlab = "Observed No. of ASVs", ylab = "Rarefied No. of ASVs")
abline(0,1)
# plot slopes
plot(S, Sslope, xlab = "Observed No. of ASVs", ylab = "Slope at rarefied sample size")
# store rarefaction curves data with 100 steps
rarecurve_data <- rarecurve(SVsControls, step = 100, sample = raremax)
# clean plot of rarefaction curves
rareplot<-map_dfr(rarecurve_data, bind_rows) %>%
bind_cols(SampleID = rownames(SVsControls),.) %>%
pivot_longer(-SampleID) %>%
drop_na() %>%
mutate(n_seqs = as.numeric(str_replace(name, "N", ""))) %>%
select(-name) %>%
left_join(metadata) %>%
# edit sampletype labels
ggplot(aes(x=n_seqs, y=value, group=SampleID, color=Group, lty=SampleType)) +
geom_line(lwd=0.6) +
scale_color_manual(values=c("lightgreen","darkgreen","lightblue","darkblue"),
labels=c("Northern Ambient", "Northern Enriched", "Southern Ambient", "Southern Enriched"))+
scale_x_continuous(labels=label_comma(), expand=c(0,0))+
scale_y_continuous(expand=c(0,0), limits=c(0,NA))+
theme(legend.key.width = unit(1.33,"cm"), legend.background = element_blank(),
legend.position="right", legend.direction = "vertical")+
labs(x = "Number of sequences", y="Number of ASVs detected", color=NULL, lty=NULL)+
guides(color=guide_legend(override.aes = list(lwd=4)))
rareplot
Summary figure
cowplot::plot_grid(reads, rareplot, align="v", axis="l", labels="AUTO", nrow=2)
ggsave("figures/PJamplicons_seqQC.png", bg="transparent", width=11, height=8)
ggsave("figures/PJamplicons_seqQC.pdf", bg="transparent", width=11, height=8)Taxa Barplots
Phylum level
# set palette
palettephy<-c("#65d6b6","#97007b","#ffbd68","#014aad","#dbacff","#ddc077",
"#f14690", "#009b75","#7a4600","grey60","grey40")
taxaAvgPhylum <-
# select aggregated phylum data
t(taxasumsclean$Phylum) %>% as.data.frame() %>%
# put rownames as SampleID
rownames_to_column("SampleID") %>%
# join with the metadata
left_join(metadata) %>%
# convert to long format
pivot_longer("d__Bacteria; Abditibacteriota":"d__Bacteria; Zixibacteria") %>%
# Clean up the phylum names: remove Bacteria; and change NA to Unknown
mutate(Phylum=str_remove(name, "d__Bacteria; "),
Phylum=recode(Phylum, "NA"="Unknown")) %>%
# Make new column PhylumOther where the least abundant phylas are grouped into "Others"
mutate(PhylumOther=forcats::fct_lump_n(f=Phylum, w=value, other_level="Others", n=10)) %>%
# calculate sums per sample per PhylumOther
group_by(SampleID, PhylumOther) %>% mutate(sumreads=sum(value)) %>%
# calculate means per bucket
group_by(SampleType, Bucket,PhylumOther) %>% mutate(bucketmeanreads=mean(sumreads)) %>%
# remove multiple values per bucket per PhylumOther
distinct(PhylumOther, Location, Treatment, SampleType, Group, Bucket, bucketmeanreads) %>%
# sum up total mean reads per phylum so we can order them
group_by(PhylumOther) %>% mutate(sumphylumOther=sum(bucketmeanreads)) %>% ungroup() %>%
# reorder PhylumOther by abundance, then move unknown and other to the end
mutate(PhylumOther=fct_reorder(PhylumOther, -sumphylumOther),
PhylumOther=fct_relevel(PhylumOther, "Unknown", after = Inf),
PhylumOther=fct_relevel(PhylumOther, "Others", after = Inf)) %>%
select(-sumphylumOther)
# start plotting
phylabarplot <- taxaAvgPhylum %>%
ggplot(aes(x=Bucket, y=bucketmeanreads, fill=PhylumOther))+
# panel per sampletype, location, treatment
facet_nested(.~SampleType+Location+Treatment, scales="free",space="free",
nest_line = element_line(color="white"),
strip = strip_nested(
text_x = elem_list_text(face = c("bold","bold","bold",NA,NA,NA,NA,NA,NA,
NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
size = c(24,24,24,20,20,20,20,20,20,
14,14,14,14,14,14,14,14,14,14,14,14))))+
# add bars
geom_col(position="fill", color="white")+
# edit the theme
theme(legend.position = "right",
legend.background = element_rect(fill=alpha("white", 0.8), color = "transparent"),
axis.text.x = element_blank(), axis.ticks.x = element_blank(),
axis.ticks.y = element_line(inherit.blank=FALSE, color="grey30"))+
# set the colors using our palette
scale_fill_manual(values=palettephy)+
# change the y-axis labels to %
scale_y_continuous(labels = scales::percent_format(), expand=c(0,0))+
# set the labels
labs(y="Percent Abundance",x=NULL,fill="Phylum")
phylabarplot
Order level
# set palette
paletteord <-c("#A6CEE3", "#1F78B4", "#B2DF8A", "#33A02C", "#FB9A99","#E31A1C", "#FDBF6F", "#FF7F00",
"#CAB2D6", "#6A3D9A", "#FFFF99", "#B15928" ,"grey60","grey40")
# select aggregated order data
t(taxasumsclean$Order) %>% as.data.frame() %>%
# put rownames as SampleID
rownames_to_column("SampleID") %>%
# join with the metadata
left_join(metadata) %>%
# convert to long format
pivot_longer("d__Bacteria; Abditibacteriota; Abditibacteria; Abditibacteriales":
"d__Bacteria; Zixibacteria; Zixibacteria; Zixibacteria") %>%
# Clean up the phylum names: remove Bacteria; and change NA to Unknown
mutate(Order=str_remove(name, "d__Bacteria; "),
Order=recode(Order, "NA; NA; NA"="Unknown")) %>%
# Make new column OrderOther where the least abundant orders are grouped into "Others"
mutate(OrderOther=forcats::fct_lump_n(f=Order, w=value, other_level="Others", n=13)) %>%
# calculate sums per sample per OrderOther
group_by(SampleID, OrderOther) %>% mutate(sumreads=sum(value)) %>%
# calculate means per bucket
group_by(SampleType, Bucket,OrderOther) %>% mutate(bucketmeanreads=mean(sumreads)) %>%
# remove multiple values per bucket per OrderOther
distinct(OrderOther, Location, Treatment, SampleType, Group, Bucket, bucketmeanreads) %>%
# sum up total mean reads per phylum so we can order them
group_by(OrderOther) %>% mutate(sumOrderOther=sum(bucketmeanreads)) %>% ungroup() %>%
# reorder OrderOther by abundance, then move unknown and other to the end
mutate(OrderOther=fct_reorder(OrderOther, -sumOrderOther),
OrderOther=fct_relevel(OrderOther, "Unknown", after = Inf),
OrderOther=fct_relevel(OrderOther, "Others", after = Inf)) %>%
# start plotting
ggplot(aes(x=Bucket, y=bucketmeanreads, fill=OrderOther))+
# panel per fish type
facet_nested(.~SampleType+Location+Treatment, scales="free",space="free",
nest_line = element_line(color="white"),
strip = strip_nested(
text_x = elem_list_text(face = c("bold","bold","bold",NA,NA,NA,NA,NA,NA,
NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
size = c(24,24,24,20,20,20,20,20,20,
14,14,14,14,14,14,14,14,14,14,14,14))))+
# add bars
geom_col(position="fill", color="white")+
# edit the theme
theme(legend.position = "bottom",
legend.background = element_rect(fill=alpha("white", 0.8), color = "transparent"),
axis.text.x = element_blank(), axis.ticks.x = element_blank(),
axis.ticks.y = element_line(inherit.blank=FALSE, color="grey30"))+
# set the colors using our palette
scale_fill_manual(values=paletteord)+
# change the y-axis labels to %
scale_y_continuous(labels = scales::percent_format(), expand=c(0,0))+
# set the labels
labs(y="Percent Abundance per Order",x=NULL,fill="Order")
# select aggregated order data
t(taxasumsclean$Order) %>% as.data.frame() %>%
# put rownames as SampleID
rownames_to_column("SampleID") %>%
# join with the metadata
left_join(metadata) %>%
# convert to long format
pivot_longer("d__Bacteria; Abditibacteriota; Abditibacteria; Abditibacteriales":
"d__Bacteria; Zixibacteria; Zixibacteria; Zixibacteria") %>%
# Clean up the phylum names: remove Bacteria; and change NA to Unknown
mutate(Order=str_remove(name, "d__Bacteria; "),
Order=recode(Order, "NA; NA; NA"="Unknown")) %>%
# Make new column OrderOther where the least abundant orders are grouped into "Others"
mutate(OrderOther=forcats::fct_lump_n(f=Order, w=value, other_level="Others", n=25)) %>%
# calculate sums per sample per OrderOther
group_by(SampleID, OrderOther) %>% mutate(sumreads=sum(value)) %>%
# calculate means per bucket
group_by(SampleType, Bucket, OrderOther) %>% mutate(bucketmeanreads=mean(sumreads)) %>%
# remove multiple values per bucket per OrderOther
distinct(OrderOther, Location, Treatment, SampleType, Group, Bucket, bucketmeanreads) %>%
# sum up total mean reads per phylum so we can order them
group_by(OrderOther) %>% mutate(sumOrderOther=sum(bucketmeanreads)) %>% ungroup() %>%
# reorder OrderOther by abundance, then move unknown and other to the end
mutate(OrderOther=fct_reorder(OrderOther, -sumOrderOther),
OrderOther=fct_relevel(OrderOther, "Unknown", after = Inf),
OrderOther=fct_relevel(OrderOther, "Others", after = Inf)) %>%
# calculate percent per OrderOther per bucket
group_by(Bucket) %>% mutate(sumBucket=sum(bucketmeanreads)) %>%
group_by(Bucket, OrderOther) %>% mutate(percOrder = bucketmeanreads/sumBucket) %>%
# Get Phylum column and order
separate_wider_delim(OrderOther, names = "Phylum", delim="; ", too_many = "drop", cols_remove=FALSE) %>%
mutate(Phylum=fct_relevel(Phylum, "Others", after = Inf),
Phylum=fct_relevel(Phylum, "Unknown", after = Inf)) %>%
# start plotting
ggplot(aes(x=Bucket, fill=percOrder, y=OrderOther))+
# panel per sampletype, location, treatment
facet_nested(Phylum~SampleType+Location+Treatment,
scales="free", space="free", nest_line = element_line(color="white"),
strip = strip_nested(
background_x = elem_list_rect(fill = alpha(c(
"orange","salmon","darkred",
"seagreen3","cornflowerblue","seagreen3","cornflowerblue","seagreen3","cornflowerblue",
"lightgreen","darkgreen","lightblue","darkblue","lightgreen","darkgreen",
"lightblue","darkblue","lightgreen","darkgreen","lightblue","darkblue"), 0.4)),
text_x = elem_list_text(face = c("bold","bold","bold",NA,NA,NA,NA,NA,NA,
NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
size = c(24,24,24,20,20,20,20,20,20,
14,14,14,14,14,14,14,14,14,14,14,14))))+
# add tiles
geom_tile()+
scale_fill_gradientn(labels = scales::percent,limits=c(0,0.4),
colours=c("white","#02818a","#016c59","#014636","black"))+
theme(axis.text.x = element_blank(), legend.position="bottom",
strip.text.x = element_text(size=16, color="black"),strip.text.y = element_blank(),
panel.border = element_rect(color="grey", fill="transparent"),
axis.ticks.y = element_line(inherit.blank=FALSE, color="grey30"),
axis.ticks.x = element_blank(), legend.background = element_blank(),
legend.key.size = unit(3, 'lines'))+
# set the labels
labs(fill="Relative Percent \nAbundance per Order",x=NULL,y=NULL)
ggsave("figures/PJamplicons_orderHeatmap.png", bg="transparent", width=16, height=9)
ggsave("figures/PJamplicons_orderHeatmap.pdf", bg="transparent", width=16, height=9)PCA with gas rates
# vector of metadata to include in the PCA
metadatavec <- metadata %>%
select(c("Bucket","Denitrification_18C","Denitrification_24C",
"NitrousOxide_18C","NitrousOxide_24C",
"August_NH4_high","August_NO3_high","August_NO2_high",
"August_NH4_low","August_NO3_low","August_NO2_low")) %>%
distinct() %>%
group_by(Bucket) %>%
# average the high/low DIN concentrations
mutate(InSitu_NH4 = mean(c(August_NH4_high, August_NH4_low), na.rm=TRUE),
InSitu_NO3 = mean(c(August_NO3_high, August_NO3_low), na.rm=TRUE),
InSitu_NO2 = mean(c(August_NO2_high, August_NO2_low), na.rm=TRUE)) %>%
select(-c("August_NH4_high","August_NO3_high","August_NO2_high",
"August_NH4_low","August_NO3_low","August_NO2_low"))
# make a table with the phylum data averaged per bucket
phytable <-
taxaAvgPhylum %>%
group_by(Bucket) %>% mutate(bucketsum=sum(bucketmeanreads)) %>%
ungroup() %>% mutate(bucketperc = bucketmeanreads/bucketsum) %>%
select(Bucket, SampleType, Group, bucketmeanreads, PhylumOther) %>% distinct() %>%
pivot_wider(names_from = PhylumOther, values_from = bucketmeanreads, values_fn = mean) %>%
# add in the metadata
left_join(metadatavec) %>%
unite("Row",Bucket, Group, SampleType) %>%
column_to_rownames("Row")
# calculate PCA
pca_result <- prcomp(phytable, scale.=TRUE, center=TRUE)
# check PCA
summary(pca_result)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.9908 1.7176 1.4588 1.32776 1.26329 1.05011 0.9703
## Proportion of Variance 0.2202 0.1639 0.1182 0.09794 0.08866 0.06126 0.0523
## Cumulative Proportion 0.2202 0.3841 0.5023 0.60024 0.68890 0.75016 0.8025
## PC8 PC9 PC10 PC11 PC12 PC13 PC14
## Standard deviation 0.91038 0.87282 0.78069 0.67201 0.4930 0.44143 0.42319
## Proportion of Variance 0.04604 0.04232 0.03386 0.02509 0.0135 0.01083 0.00995
## Cumulative Proportion 0.84851 0.89083 0.92469 0.94978 0.9633 0.97411 0.98406
## PC15 PC16 PC17 PC18
## Standard deviation 0.38558 0.26910 0.23066 0.11269
## Proportion of Variance 0.00826 0.00402 0.00296 0.00071
## Cumulative Proportion 0.99232 0.99634 0.99929 1.00000# get metadata averaged per bucket
metadatabucket <- metadata %>%
unite("Row",Bucket, Group, SampleType, remove=FALSE) %>%
select(-c(SampleName, SampleID, OysterNumber)) %>%
distinct() %>%
drop_na(SampleType)
library(ggfortify)
# plot PCA
aplot<-autoplot(pca_result,size=6, data=metadatabucket,
fill="TypeGroup",colour="TypeGroup",
shape="SampleType",
loadings.label.repel=TRUE,loadings.label.size=6,
loadings.label.colour=c("grey30","grey30","grey30","grey30","grey30",
"grey30","grey30","grey30","grey30","grey30","grey30",
"red","maroon","red","maroon","purple","purple","purple"),
loadings.colour=c("grey30","grey30","grey30","grey30","grey30",
"grey30","grey30","grey30","grey30","grey30","grey30",
"red","maroon","red","maroon","purple","purple","purple"),
loadings=TRUE,loadings.label=TRUE,
frame=TRUE,frame.colour = "TypeGroup")
PCAplot <- aplot + scale_fill_manual(values=c("lightgreen","darkgreen","lightblue","darkblue",
"lightgreen","darkgreen","lightblue","darkblue",
"lightgreen","darkgreen","lightblue","darkblue")) +
scale_color_manual(values=c("lightgreen","darkgreen","lightblue","darkblue",
"lightgreen","darkgreen","lightblue","darkblue",
"lightgreen","darkgreen","lightblue","darkblue")) +
theme(legend.position="none")
PCAlegend <- get_legend(ggplot(metadatabucket, aes(x=Bucket, y=Bucket, shape=SampleType, color=Group))+
geom_point(size=6)+
scale_color_manual(values=c("lightgreen","darkgreen","lightblue","darkblue"),
labels=c("Northern Ambient", "Northern Enriched",
"Southern Ambient", "Southern Enriched"))+
theme(legend.position = "right")+
guides(color=guide_legend(override.aes = list(size=8, shape=15))))
PCAplot+PCAlegend+plot_layout(widths = c(3,1))
ggsave("figures/PJamplicons_PCA.png", bg="transparent", width=10, height=9)
ggsave("figures/PJamplicons_PCA.pdf", bg="transparent", width=8, height=10)Venn Diagrams
SampleType
At the ASV level for gut/inner/outer
# list of gut ASVs
gutset <- dataASVsClean %>% ungroup() %>% filter(SampleType=="Gut") %>%
filter(value!=0) %>% distinct(name) %>% .$name
# list of outer ASVs
outerset <- dataASVsClean %>% ungroup() %>% filter(SampleType=="Outer Shell") %>%
filter(value!=0) %>% distinct(name) %>% .$name
# list of inner ASVs
innerset <- dataASVsClean %>% ungroup() %>% filter(SampleType=="Inner Shell") %>%
filter(value!=0) %>% distinct(name) %>% .$name
setsType <- list("Gut" = gutset, "Inner Shell" = innerset, "Outer Shell" = outerset)
# plot venn diagram
vennType <- ggVennDiagram::ggVennDiagram(setsType,set_size = 8,
label_alpha = 0,label_size = 6)+
scale_color_manual(values = c("orange", "tomato", "darkred"))+
scale_fill_gradient(low="white",high="grey50")+
coord_sf(clip = 'off')+theme(legend.position = "none")
vennType
library(UpSetR)
mat <- make_comb_mat(setsType)
upsetplotType <- UpSet(mat,
pt_size = unit(.5, "cm"),lwd=2.5,
comb_col = c("black","grey50","grey50","grey50","orange", "darkred", "tomato"),
left_annotation = upset_left_annotation(mat, bar_width=0.7,
axis_param = list(side = "bottom",labels_rot = 0),
# annotation_name_side = "top",
gp = gpar(fill = c("orange", "darkred", "tomato")),
width = unit(4, "cm")),
row_names_side = "left",
top_annotation = upset_top_annotation(mat,bar_width = 0.9, height = unit(6, "cm")))
upsetplotType
Site and enrichment
Northern_Ambientset <- dataASVsClean %>% ungroup() %>% filter(Group=="Northern_Ambient") %>%
filter(value!=0) %>% distinct(name) %>% .$name
Northern_Enrichedset <- dataASVsClean %>% ungroup() %>% filter(Group=="Northern_Enriched") %>%
filter(value!=0) %>% distinct(name) %>% .$name
Southern_Ambientset <- dataASVsClean %>% ungroup() %>% filter(Group=="Southern_Ambient") %>%
filter(value!=0) %>% distinct(name) %>% .$name
Southern_Enrichedset <- dataASVsClean %>% ungroup() %>% filter(Group=="Southern_Enriched") %>%
filter(value!=0) %>% distinct(name) %>% .$name
setsSiteTreat <- list("Northern Ambient" = Northern_Ambientset, "Northern Enriched" = Northern_Enrichedset,
"Southern Ambient" = Southern_Ambientset, "Southern Enriched" = Southern_Enrichedset)
# plot venn diagram
vennSiteTreat <- ggVennDiagram::ggVennDiagram(setsSiteTreat,set_size = 8,
label_alpha = 0,label_size = 6)+
scale_color_manual(values = c("lightgreen","darkgreen","lightblue","darkblue"))+
scale_fill_gradient(low="white",high="grey50")+
coord_sf(clip = 'off')+theme(legend.position = "none")
vennSiteTreat
matST <- make_comb_mat(setsSiteTreat)
upsetSiteTreat <- UpSet(matST,
pt_size = unit(.5, "cm"),lwd=2.5,
comb_col = c("black","grey50","grey50","grey50","grey50",
"cornflowerblue","grey50","grey50","grey50","grey50",
"seagreen3","lightgreen","darkgreen","lightblue","darkblue"),
left_annotation = upset_left_annotation(matST, bar_width=0.7,
axis_param = list(side = "bottom",labels_rot = 0),
# annotation_name_side = "top",
gp = gpar(fill = c("lightgreen","darkgreen","lightblue","darkblue")),
width = unit(4, "cm")),
row_names_side = "left",
top_annotation = upset_top_annotation(matST,bar_width = 0.9, height = unit(6, "cm")))
upsetSiteTreat
Summary figure
cowplot::plot_grid(ggplotify::as.grob(upsetplotType),
ggplotify::as.grob(upsetSiteTreat), nrow=2, labels="AUTO")
ggsave("figures/PJamplicons_upsetPlots.png", bg="transparent", width=7, height=8)
ggsave("figures/PJamplicons_upsetPlots.pdf", bg="transparent", width=8, height=8)Alpha Diversity
metasamples <- metadata %>% drop_na(Denitrification_18C)
# calculate chao diversity
chao <- estimateR(SVsControls) %>% t() %>% as.data.frame()
# calculate shannon
diversityshannon<-diversity(SVsControls, index="shannon")
# calculate simpson
diversitysimpsons<-diversity(SVsControls, index="simpson")
# add into metadata variable
metasamples$Chao <- chao$S.chao1
metasamples$Simpsons<-diversitysimpsons
metasamples$Shannon<-diversityshannon
metasamples %>%
pivot_longer(Chao:Shannon, names_to="DiversityIndex", values_to="Value") %>%
ggplot(aes(x=Group,y=Value, fill=Group, shape=Group))+
geom_jitter(width=0.15, size=3, alpha=0.8)+
geom_boxplot(alpha=0.8)+
geom_hline(yintercept = 0)+
stat_compare_means(label.y=0.3)+
facet_grid(DiversityIndex~SampleType, scales="free")+
labs(x=NULL, y="Index of Diversity")+
scale_fill_manual(values=c("lightgreen","darkgreen","lightblue","darkblue"))+
scale_shape_manual(values=c(22,23,24,25))+
scale_y_continuous(expand = expansion(mult = c(0, .1)), limits=c(0,NA))+
theme(legend.position = "bottom",
axis.text.x = element_blank())
Beta Diversity
# ellipse function
veganCovEllipse<-function (cov, center = c(0, 0), scale = 1, npoints = 100)
{theta <- (0:npoints) * 2 * pi/npoints
Circle <- cbind(cos(theta), sin(theta))
t(center + scale * t(Circle %*% chol(cov)))}All samples together
# metadata without controls
metadataoyster <- metadata %>% drop_na(Denitrification_18C) %>%
mutate(Group=recode(Group,
"Northern_Ambient"="Northern Ambient",
"Northern_Enriched"="Northern Enriched",
"Southern_Ambient"="Southern Ambient",
"Southern_Enriched"="Southern Enriched"))
# calculate beta-diversity
sol<-metaMDS(SVsControls,distance = "bray", k = 2, trymax = 200)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1479886
## Run 1 stress 0.1747572
## Run 2 stress 0.1714223
## Run 3 stress 0.1522169
## Run 4 stress 0.1540459
## Run 5 stress 0.1540459
## Run 6 stress 0.1484422
## ... Procrustes: rmse 0.007908788 max resid 0.07878964
## Run 7 stress 0.1793566
## Run 8 stress 0.1481364
## ... Procrustes: rmse 0.004402954 max resid 0.04220882
## Run 9 stress 0.1572643
## Run 10 stress 0.1835301
## Run 11 stress 0.1870935
## Run 12 stress 0.1572643
## Run 13 stress 0.1479886
## ... New best solution
## ... Procrustes: rmse 7.528299e-05 max resid 0.0006525791
## ... Similar to previous best
## Run 14 stress 0.1481365
## ... Procrustes: rmse 0.004401011 max resid 0.0423366
## Run 15 stress 0.1647826
## Run 16 stress 0.1647826
## Run 17 stress 0.176382
## Run 18 stress 0.1871848
## Run 19 stress 0.1899911
## Run 20 stress 0.1782796
## *** Best solution repeated 1 times# plot the solution
sol$stress; stressplot(sol)## [1] 0.1479886# make clean data with solution
NMDS=data.frame(x=sol$point[,1],y=sol$point[,2],
SampleType=as.factor(metadataoyster$SampleType),
Group=as.factor(metadataoyster$Group))by sample type
# make ellipse paths with solution
plot.new()
ordType<-ordiellipse(sol, NMDS$SampleType, display = "sites", kind ="sd", conf = 0.95, label = T)dev.off()## null device
## 1# store ellipse paths
df_ellType <- data.frame()
for(g in levels(NMDS$SampleType)){
if(g!="" && (g %in% names(ordType))){
df_ellType <- rbind(df_ellType, cbind(as.data.frame(with(NMDS[NMDS$SampleType==g,],
veganCovEllipse(ordType[[g]]$cov,ordType[[g]]$center,ordType[[g]]$scale))),SampleType=g))}}
head(df_ellType)## NMDS1 NMDS2 SampleType
## 1 -0.2272305 0.4340363 Gut
## 2 -0.2285191 0.5001351 Gut
## 3 -0.2323800 0.5645280 Gut
## 4 -0.2387979 0.6269609 Gut
## 5 -0.2477475 0.6871874 Gut
## 6 -0.2591935 0.7449698 Gut# determine center of ellipses
NMDS.mean.type=aggregate(NMDS[,1:2],list(group=NMDS$SampleType),mean)
adonis2(SVsControls~SampleType, data=metadataoyster, by=NULL,method="bray", k=2)## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = SVsControls ~ SampleType, data = metadataoyster, method = "bray", by = NULL, k = 2)
## Df SumOfSqs R2 F Pr(>F)
## Model 2 7.922 0.24702 17.223 0.001 ***
## Residual 105 24.149 0.75298
## Total 107 32.071 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1betadivType <- ggplot(data=NMDS,aes(x,y,colour=SampleType, fill=SampleType))+
geom_path(data=df_ellType, aes(x=NMDS1, y=NMDS2, lty=SampleType), size=1) +
geom_point(size=4, alpha=0.8, aes(shape=Group))+
scale_shape_manual(values=c(22,23,24,25))+
annotate("text",x=NMDS.mean.type$x,y=NMDS.mean.type$y,
label=NMDS.mean.type$group, size=5, color="gray40") +
scale_fill_manual(values=c("orange", "tomato", "darkred"))+
scale_colour_manual(values=c("orange", "tomato", "darkred"))+
labs(x=NULL, y=NULL)+
ggtitle("Sample Type")+
theme(legend.position="right")
betadivType
by location and enrichment
# make ellipse paths with solution
plot.new()
ordGroup<-ordiellipse(sol, NMDS$Group, display = "sites", kind ="sd", conf = 0.95, label = T)dev.off()## null device
## 1# store ellipse paths
df_ellGroup <- data.frame()
for(g in levels(NMDS$Group)){
if(g!="" && (g %in% names(ordGroup))){
df_ellGroup <- rbind(df_ellGroup, cbind(as.data.frame(with(NMDS[NMDS$Group==g,],
veganCovEllipse(ordGroup[[g]]$cov,ordGroup[[g]]$center,ordGroup[[g]]$scale))),Group=g))}}
head(df_ellGroup)## NMDS1 NMDS2 Group
## 1 1.776553 0.3835542 Northern Ambient
## 2 1.773278 0.4050816 Northern Ambient
## 3 1.763466 0.4265556 Northern Ambient
## 4 1.747155 0.4478915 Northern Ambient
## 5 1.724410 0.4690052 Northern Ambient
## 6 1.695321 0.4898132 Northern Ambient# determine center of ellipses
NMDS.mean.group=aggregate(NMDS[,1:2],list(group=NMDS$Group),mean)
adonis2(SVsControls~Group, data=metadataoyster, by=NULL,method="bray", k=2)## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = SVsControls ~ Group, data = metadataoyster, method = "bray", by = NULL, k = 2)
## Df SumOfSqs R2 F Pr(>F)
## Model 3 3.818 0.11905 4.6849 0.001 ***
## Residual 104 28.252 0.88095
## Total 107 32.071 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1betadivGroup <- ggplot(data=NMDS,aes(x,y,colour=Group, fill=Group))+
geom_path(data=df_ellGroup, aes(x=NMDS1, y=NMDS2, lty=Group), size=1) +
geom_point(size=4, alpha=0.8, aes(shape=SampleType))+
scale_shape_manual(values=c(21,24,22))+
annotate("text",x=NMDS.mean.group$x,y=NMDS.mean.group$y,
label=NMDS.mean.group$group, size=5, color="gray40") +
scale_color_manual(values=c("lightgreen","darkgreen","lightblue","darkblue"))+
scale_fill_manual(values=c("lightgreen","darkgreen","lightblue","darkblue"))+
labs(x=NULL, y=NULL)+
ggtitle("Site and Treatment")+
theme(legend.position="right")
betadivGroup
Separate by sample type
Gut
# gut sample metadata
metadataoyster_gut <- metadataoyster %>% filter(SampleType=="Gut")
# gut sample ASVs
SVs_gut <- as.data.frame(SVsControls) %>% rownames_to_column("SampleID") %>%
filter(SampleID %in% metadataoyster_gut$SampleID) %>%
column_to_rownames("SampleID") %>% as.matrix()
# calculate beta-diversity
sol_gut<-metaMDS(SVs_gut,distance = "bray", k = 2, trymax = 200)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.152747
## Run 1 stress 0.1768887
## Run 2 stress 0.1527469
## ... New best solution
## ... Procrustes: rmse 9.108185e-05 max resid 0.000340311
## ... Similar to previous best
## Run 3 stress 0.1536736
## Run 4 stress 0.184936
## Run 5 stress 0.1527469
## ... Procrustes: rmse 6.673442e-05 max resid 0.0002724732
## ... Similar to previous best
## Run 6 stress 0.1791047
## Run 7 stress 0.1793466
## Run 8 stress 0.1527469
## ... Procrustes: rmse 9.896829e-05 max resid 0.0005100485
## ... Similar to previous best
## Run 9 stress 0.152747
## ... Procrustes: rmse 0.0001406771 max resid 0.0006561834
## ... Similar to previous best
## Run 10 stress 0.1741098
## Run 11 stress 0.1527469
## ... New best solution
## ... Procrustes: rmse 5.934536e-06 max resid 2.12092e-05
## ... Similar to previous best
## Run 12 stress 0.1527469
## ... Procrustes: rmse 7.086661e-05 max resid 0.0002953921
## ... Similar to previous best
## Run 13 stress 0.1763947
## Run 14 stress 0.1743605
## Run 15 stress 0.1781819
## Run 16 stress 0.1527469
## ... Procrustes: rmse 0.0001405581 max resid 0.000686137
## ... Similar to previous best
## Run 17 stress 0.1795616
## Run 18 stress 0.1886936
## Run 19 stress 0.1527469
## ... Procrustes: rmse 3.489591e-05 max resid 0.0001323789
## ... Similar to previous best
## Run 20 stress 0.1743605
## *** Best solution repeated 4 times# plot the solution
sol_gut$stress; stressplot(sol_gut)## [1] 0.1527469# make clean data with solution
NMDS_gut=data.frame(x=sol_gut$point[,1],y=sol_gut$point[,2],
Group=as.factor(metadataoyster_gut$Group))
# make ellipse paths with solution
plot.new()
ord_gut<-ordiellipse(sol_gut, NMDS_gut$Group, display = "sites", kind ="sd", conf = 0.95, label = T)dev.off()## null device
## 1# store ellipse paths
df_ell_gut <- data.frame()
for(g in levels(NMDS$Group)){
if(g!="" && (g %in% names(ord_gut))){
df_ell_gut <- rbind(df_ell_gut, cbind(as.data.frame(with(NMDS_gut[NMDS_gut$Group==g,],
veganCovEllipse(ord_gut[[g]]$cov,ord_gut[[g]]$center,ord_gut[[g]]$scale))),Group=g))}}
head(df_ell_gut)## NMDS1 NMDS2 Group
## 1 0.4521084 -0.2410993 Northern Ambient
## 2 0.4516529 -0.2273899 Northern Ambient
## 3 0.4502880 -0.2129048 Northern Ambient
## 4 0.4480193 -0.1977010 Northern Ambient
## 5 0.4448556 -0.1818387 Northern Ambient
## 6 0.4408094 -0.1653803 Northern Ambientadonis2(SVs_gut~Group, data=metadataoyster_gut, by=NULL,method="bray", k=2)## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = SVs_gut ~ Group, data = metadataoyster_gut, method = "bray", by = NULL, k = 2)
## Df SumOfSqs R2 F Pr(>F)
## Model 3 1.8145 0.22457 3.0892 0.001 ***
## Residual 32 6.2654 0.77543
## Total 35 8.0799 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1betadiv_gut <- ggplot(data=NMDS_gut,aes(x,y,colour=Group, fill=Group))+
geom_path(data=df_ell_gut, aes(x=NMDS1, y=NMDS2, lty=Group), size=1) +
geom_point(size=4, alpha=0.8, aes(shape=Group))+
scale_shape_manual(values=c(21,22,21,22))+
scale_color_manual(values=c("lightgreen","darkgreen","lightblue","darkblue"))+
scale_fill_manual(values=c("lightgreen","darkgreen","lightblue","darkblue"))+
labs(x=NULL, y=NULL, title="Gut")+
theme(legend.position="right")
betadiv_gut
Inner shell
# inner sample metadata
metadataoyster_inner <- metadataoyster %>% filter(SampleType=="Inner Shell")
# inner sample ASVs
SVs_inner <- as.data.frame(SVsControls) %>% rownames_to_column("SampleID") %>%
filter(SampleID %in% metadataoyster_inner$SampleID) %>%
column_to_rownames("SampleID") %>% as.matrix()
# calculate beta-diversity
sol_inner<-metaMDS(SVs_inner,distance = "bray", k = 2, trymax = 200)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1403371
## Run 1 stress 0.1504228
## Run 2 stress 0.1647692
## Run 3 stress 0.1403437
## ... Procrustes: rmse 0.001293525 max resid 0.005918181
## ... Similar to previous best
## Run 4 stress 0.1504232
## Run 5 stress 0.140858
## Run 6 stress 0.1415385
## Run 7 stress 0.1636471
## Run 8 stress 0.1403436
## ... Procrustes: rmse 0.001297858 max resid 0.005916375
## ... Similar to previous best
## Run 9 stress 0.1412514
## Run 10 stress 0.1412067
## Run 11 stress 0.1412515
## Run 12 stress 0.1403369
## ... New best solution
## ... Procrustes: rmse 0.0001244175 max resid 0.0003872171
## ... Similar to previous best
## Run 13 stress 0.140344
## ... Procrustes: rmse 0.001441538 max resid 0.006084874
## ... Similar to previous best
## Run 14 stress 0.3871885
## Run 15 stress 0.1414897
## Run 16 stress 0.1678357
## Run 17 stress 0.1414898
## Run 18 stress 0.1412514
## Run 19 stress 0.1412972
## Run 20 stress 0.1504232
## *** Best solution repeated 2 times# plot the solution
sol_inner$stress; stressplot(sol_inner)## [1] 0.1403369# make clean data with solution
NMDS_inner=data.frame(x=sol_inner$point[,1],y=sol_inner$point[,2],
Group=as.factor(metadataoyster_inner$Group))
# make ellipse paths with solution
plot.new()
ord_inner<-ordiellipse(sol_inner, NMDS_inner$Group, display = "sites", kind ="sd", conf = 0.95, label = T)dev.off()## null device
## 1# store ellipse paths
df_ell_inner <- data.frame()
for(g in levels(NMDS$Group)){
if(g!="" && (g %in% names(ord_inner))){
df_ell_inner <- rbind(df_ell_inner, cbind(as.data.frame(with(NMDS_inner[NMDS_inner$Group==g,],
veganCovEllipse(ord_inner[[g]]$cov,ord_inner[[g]]$center,ord_inner[[g]]$scale))),Group=g))}}
head(df_ell_inner)## NMDS1 NMDS2 Group
## 1 0.9328061 0.009293697 Northern Ambient
## 2 0.9317115 0.023472774 Northern Ambient
## 3 0.9284321 0.038304847 Northern Ambient
## 4 0.9229809 0.053731382 Northern Ambient
## 5 0.9153793 0.069691498 Northern Ambient
## 6 0.9056573 0.086122206 Northern Ambientadonis2(SVs_inner~Group, data=metadataoyster_inner, by=NULL,method="bray", k=2)## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = SVs_inner ~ Group, data = metadataoyster_inner, method = "bray", by = NULL, k = 2)
## Df SumOfSqs R2 F Pr(>F)
## Model 3 1.6891 0.18767 2.4642 0.001 ***
## Residual 32 7.3113 0.81233
## Total 35 9.0003 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1betadiv_inner <- ggplot(data=NMDS_inner,aes(x,y,colour=Group, fill=Group))+
geom_path(data=df_ell_inner, aes(x=NMDS1, y=NMDS2, lty=Group), size=1) +
geom_point(size=4, alpha=0.8, aes(shape=Group))+
scale_shape_manual(values=c(21,22,21,22))+
scale_color_manual(values=c("lightgreen","darkgreen","lightblue","darkblue"))+
scale_fill_manual(values=c("lightgreen","darkgreen","lightblue","darkblue"))+
labs(x=NULL, y=NULL, title="Inner Shell")+
theme(legend.position="right")
betadiv_inner
Outer shell
# outer sample metadata
metadataoyster_outer <- metadataoyster %>% filter(SampleType=="Outer Shell")
# outer sample ASVs
SVs_outer <- as.data.frame(SVsControls) %>% rownames_to_column("SampleID") %>%
filter(SampleID %in% metadataoyster_outer$SampleID) %>%
column_to_rownames("SampleID") %>% as.matrix()
# calculate beta-diversity
sol_outer<-metaMDS(SVs_outer,distance = "bray", k = 2, trymax = 200)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.07049937
## Run 1 stress 0.07049946
## ... Procrustes: rmse 0.001026726 max resid 0.005484579
## ... Similar to previous best
## Run 2 stress 0.07049883
## ... New best solution
## ... Procrustes: rmse 0.0007120391 max resid 0.003792605
## ... Similar to previous best
## Run 3 stress 0.07049928
## ... Procrustes: rmse 0.0002503437 max resid 0.001341534
## ... Similar to previous best
## Run 4 stress 0.0704993
## ... Procrustes: rmse 0.0002541802 max resid 0.001362108
## ... Similar to previous best
## Run 5 stress 0.0704994
## ... Procrustes: rmse 0.0002836352 max resid 0.001518798
## ... Similar to previous best
## Run 6 stress 0.07049927
## ... Procrustes: rmse 0.0002424862 max resid 0.001298387
## ... Similar to previous best
## Run 7 stress 0.07049918
## ... Procrustes: rmse 0.0002088187 max resid 0.001117844
## ... Similar to previous best
## Run 8 stress 0.07049898
## ... Procrustes: rmse 0.000107109 max resid 0.0005708784
## ... Similar to previous best
## Run 9 stress 0.07049906
## ... Procrustes: rmse 0.0005737213 max resid 0.00305365
## ... Similar to previous best
## Run 10 stress 0.07049907
## ... Procrustes: rmse 0.0001500068 max resid 0.0008007891
## ... Similar to previous best
## Run 11 stress 0.07049965
## ... Procrustes: rmse 0.0008099288 max resid 0.004317483
## ... Similar to previous best
## Run 12 stress 0.0704989
## ... Procrustes: rmse 5.147471e-05 max resid 0.0002711369
## ... Similar to previous best
## Run 13 stress 0.07049882
## ... New best solution
## ... Procrustes: rmse 1.82649e-05 max resid 9.414951e-05
## ... Similar to previous best
## Run 14 stress 0.07049939
## ... Procrustes: rmse 0.0003014988 max resid 0.001616179
## ... Similar to previous best
## Run 15 stress 0.07049936
## ... Procrustes: rmse 0.0003010603 max resid 0.001613476
## ... Similar to previous best
## Run 16 stress 0.07049935
## ... Procrustes: rmse 0.0002877997 max resid 0.001540151
## ... Similar to previous best
## Run 17 stress 0.07049921
## ... Procrustes: rmse 0.0002416272 max resid 0.001293612
## ... Similar to previous best
## Run 18 stress 0.07049895
## ... Procrustes: rmse 8.962672e-05 max resid 0.0004669461
## ... Similar to previous best
## Run 19 stress 0.07049934
## ... Procrustes: rmse 0.0002929964 max resid 0.001570201
## ... Similar to previous best
## Run 20 stress 0.07049881
## ... New best solution
## ... Procrustes: rmse 1.142461e-05 max resid 4.530801e-05
## ... Similar to previous best
## *** Best solution repeated 1 times# plot the solution
sol_outer$stress; stressplot(sol_outer)## [1] 0.07049881# make clean data with solution
NMDS_outer=data.frame(x=sol_outer$point[,1],y=sol_outer$point[,2],
Group=as.factor(metadataoyster_outer$Group))
# make ellipse paths with solution
plot.new()
ord_outer<-ordiellipse(sol_outer, NMDS_outer$Group, display = "sites", kind ="sd", conf = 0.95, label = T)dev.off()## null device
## 1# store ellipse paths
df_ell_outer <- data.frame()
for(g in levels(NMDS$Group)){
if(g!="" && (g %in% names(ord_outer))){
df_ell_outer <- rbind(df_ell_outer, cbind(as.data.frame(with(NMDS_outer[NMDS_outer$Group==g,],
veganCovEllipse(ord_outer[[g]]$cov,ord_outer[[g]]$center,ord_outer[[g]]$scale))),Group=g))}}
head(df_ell_outer)## NMDS1 NMDS2 Group
## 1 0.8263029 -0.075594958 Northern Ambient
## 2 0.8257595 -0.054854028 Northern Ambient
## 3 0.8241311 -0.034565863 Northern Ambient
## 4 0.8214244 -0.014810529 Northern Ambient
## 5 0.8176500 0.004334006 Northern Ambient
## 6 0.8128227 0.022792190 Northern Ambientadonis2(SVs_outer~Group, data=metadataoyster_outer, by=NULL,method="bray", k=2)## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = SVs_outer ~ Group, data = metadataoyster_outer, method = "bray", by = NULL, k = 2)
## Df SumOfSqs R2 F Pr(>F)
## Model 3 3.5839 0.50704 10.971 0.001 ***
## Residual 32 3.4844 0.49296
## Total 35 7.0683 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1betadiv_outer <- ggplot(data=NMDS_outer,aes(x,y,colour=Group, fill=Group))+
geom_path(data=df_ell_outer, aes(x=NMDS1, y=NMDS2, lty=Group), size=1) +
geom_point(size=4, alpha=0.8, aes(shape=Group))+
scale_shape_manual(values=c(21,22,21,22))+
scale_color_manual(values=c("lightgreen","darkgreen","lightblue","darkblue"))+
scale_fill_manual(values=c("lightgreen","darkgreen","lightblue","darkblue"))+
labs(x=NULL, y=NULL, title="Outer Shell")+
theme(legend.position="right")
betadiv_outer
Summary figure
betadiv_gut+betadiv_inner+betadiv_outer + guide_area() + plot_layout(guides="collect") +
plot_annotation(tag_levels = "A") &
theme(legend.position="right", plot.tag = element_text(size=20))
Beta Diversity - per bucket
# metadata without controls
metadataoyster <- metadatabucket %>% drop_na(Denitrification_18C) %>%
unite("BucketType", Bucket, SampleType, remove=FALSE) %>%
mutate(Group=recode(Group,
"Northern_Ambient"="Northern Ambient",
"Northern_Enriched"="Northern Enriched",
"Southern_Ambient"="Southern Ambient",
"Southern_Enriched"="Southern Enriched"))
SVsControlsAvgFormat <- as.data.frame(SVsControlsAvg) %>% rownames_to_column("BucketType") %>%
# reorder based on metadata
slice(match(metadataoyster$BucketType, BucketType)) %>%
# convert to matrix
column_to_rownames("BucketType") %>% as.matrix()
# calculate beta-diversity
sol<-metaMDS(SVsControlsAvgFormat,distance = "bray", k = 2, trymax = 200)## Run 0 stress 0.1070776
## Run 1 stress 0.1198059
## Run 2 stress 0.1134796
## Run 3 stress 0.1069693
## ... New best solution
## ... Procrustes: rmse 0.003045396 max resid 0.01202749
## Run 4 stress 0.1198059
## Run 5 stress 0.1516368
## Run 6 stress 0.119884
## Run 7 stress 0.119884
## Run 8 stress 0.1198084
## Run 9 stress 0.1070776
## ... Procrustes: rmse 0.003049885 max resid 0.01207713
## Run 10 stress 0.1069693
## ... New best solution
## ... Procrustes: rmse 1.338803e-06 max resid 3.884138e-06
## ... Similar to previous best
## Run 11 stress 0.1614398
## Run 12 stress 0.1069693
## ... New best solution
## ... Procrustes: rmse 1.997572e-06 max resid 6.585589e-06
## ... Similar to previous best
## Run 13 stress 0.1450137
## Run 14 stress 0.10697
## ... Procrustes: rmse 0.001434899 max resid 0.006328054
## ... Similar to previous best
## Run 15 stress 0.1134796
## Run 16 stress 0.1069693
## ... Procrustes: rmse 2.984309e-06 max resid 1.1196e-05
## ... Similar to previous best
## Run 17 stress 0.1069693
## ... Procrustes: rmse 1.713126e-06 max resid 5.174092e-06
## ... Similar to previous best
## Run 18 stress 0.1291109
## Run 19 stress 0.10697
## ... Procrustes: rmse 0.001434709 max resid 0.00632356
## ... Similar to previous best
## Run 20 stress 0.1069693
## ... Procrustes: rmse 5.144497e-06 max resid 1.666774e-05
## ... Similar to previous best
## *** Best solution repeated 6 times# plot the solution
sol$stress; stressplot(sol)## [1] 0.1069693# make clean data with solution
NMDS=data.frame(x=sol$point[,1],y=sol$point[,2],
SampleType=as.factor(metadataoyster$SampleType),
Group=as.factor(metadataoyster$Group))by sample type
# make ellipse paths with solution
plot.new()
ordType<-ordiellipse(sol, NMDS$SampleType, display = "sites", kind ="sd", conf = 0.95, label = T)dev.off()## null device
## 1# store ellipse paths
df_ellType <- data.frame()
for(g in levels(NMDS$SampleType)){
if(g!="" && (g %in% names(ordType))){
df_ellType <- rbind(df_ellType, cbind(as.data.frame(with(NMDS[NMDS$SampleType==g,],
veganCovEllipse(ordType[[g]]$cov,ordType[[g]]$center,ordType[[g]]$scale))),SampleType=g))}}
head(df_ellType)## NMDS1 NMDS2 SampleType
## 1 -0.4599791 0.1006204 Gut
## 2 -0.4607695 0.1477173 Gut
## 3 -0.4631377 0.1943754 Gut
## 4 -0.4670743 0.2404106 Gut
## 5 -0.4725638 0.2856413 Gut
## 6 -0.4795846 0.3298890 Gut# determine center of ellipses
NMDS.mean.type=aggregate(NMDS[,1:2],list(group=NMDS$SampleType),mean)
adonis2(SVsControlsAvgFormat~SampleType, data=metadataoyster, by=NULL,method="bray", k=2)## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = SVsControlsAvgFormat ~ SampleType, data = metadataoyster, method = "bray", by = NULL, k = 2)
## Df SumOfSqs R2 F Pr(>F)
## Model 2 3.4487 0.39024 10.56 0.001 ***
## Residual 33 5.3887 0.60976
## Total 35 8.8375 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1betadivType <- ggplot(data=NMDS,aes(x,y,colour=SampleType, fill=SampleType))+
geom_path(data=df_ellType, aes(x=NMDS1, y=NMDS2, lty=SampleType), size=1) +
geom_point(size=4, alpha=0.8, aes(shape=Group))+
scale_shape_manual(values=c(22,23,24,25))+
annotate("text",x=NMDS.mean.type$x,y=NMDS.mean.type$y,
label=NMDS.mean.type$group, size=5,
color=c("orange", "tomato", "darkred")) +
scale_fill_manual(values=c("orange", "tomato", "darkred"))+
scale_colour_manual(values=c("orange", "tomato", "darkred"))+
labs(x=NULL, y=NULL)+
theme(legend.position="right")
betadivType
by location and enrichment
# make ellipse paths with solution
plot.new()
ordGroup<-ordiellipse(sol, NMDS$Group, display = "sites", kind ="sd", conf = 0.95, label = T)dev.off()## null device
## 1# store ellipse paths
df_ellGroup <- data.frame()
for(g in levels(NMDS$Group)){
if(g!="" && (g %in% names(ordGroup))){
df_ellGroup <- rbind(df_ellGroup, cbind(as.data.frame(with(NMDS[NMDS$Group==g,],
veganCovEllipse(ordGroup[[g]]$cov,ordGroup[[g]]$center,ordGroup[[g]]$scale))),Group=g))}}
head(df_ellGroup)## NMDS1 NMDS2 Group
## 1 1.891943 0.4774974 Northern Ambient
## 2 1.888455 0.4856031 Northern Ambient
## 3 1.878004 0.4930745 Northern Ambient
## 4 1.860633 0.4998819 Northern Ambient
## 5 1.836408 0.5059985 Northern Ambient
## 6 1.805426 0.5114001 Northern Ambient# determine center of ellipses
NMDS.mean.group=aggregate(NMDS[,1:2],list(group=NMDS$Group),mean)
adonis2(SVsControlsAvgFormat~Group, data=metadataoyster, by=NULL,method="bray", k=2)## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = SVsControlsAvgFormat ~ Group, data = metadataoyster, method = "bray", by = NULL, k = 2)
## Df SumOfSqs R2 F Pr(>F)
## Model 3 1.5792 0.1787 2.3209 0.002 **
## Residual 32 7.2582 0.8213
## Total 35 8.8375 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1betadivGroup <- ggplot(data=NMDS,aes(x,y,color=Group, fill=Group))+
geom_path(data=df_ellGroup, aes(x=NMDS1, y=NMDS2, lty=Group), size=1) +
geom_point(size=4, alpha=0.8, aes(shape=SampleType))+
scale_shape_manual(values=c(21,24,22))+
annotate("text",x=NMDS.mean.group$x,y=NMDS.mean.group$y,
label=NMDS.mean.group$group, size=5,
color=c("lightgreen","darkgreen","lightblue","darkblue")) +
scale_color_manual(values=c("lightgreen","darkgreen","lightblue","darkblue"))+
scale_fill_manual(values=c("lightgreen","darkgreen","lightblue","darkblue"))+
labs(x=NULL, y=NULL)+
theme(legend.position="right")
betadivGroup
Summary figure
phylabarplot /
(betadivType+betadivGroup) +
plot_annotation(tag_levels = "A") &
theme(legend.position="right", plot.tag = element_text(size=20))
ggsave("figures/DiversitySummary.png", bg="transparent", width=14, height=10)
ggsave("figures/DiversitySummary.pdf", bg="transparent", width=14, height=10)sessionInfo()## R version 4.1.3 (2022-03-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur/Monterey 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] ggfortify_0.4.15 ComplexHeatmap_2.10.0 ggh4x_0.2.3.9000
## [4] gplots_3.1.3 UpSetR_1.4.0 patchwork_1.1.2
## [7] gt_0.8.0 vegan_2.6-4 lattice_0.20-45
## [10] permute_0.9-7 qiime2R_0.99.6 ggpubr_0.6.0
## [13] scales_1.2.1 lubridate_1.9.2 forcats_1.0.0
## [16] stringr_1.5.0 dplyr_1.1.0 purrr_1.0.1
## [19] readr_2.1.4 tidyr_1.3.0 tibble_3.1.8
## [22] ggplot2_3.4.1 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] backports_1.4.1 circlize_0.4.15 Hmisc_4.8-0
## [4] systemfonts_1.0.4 plyr_1.8.8 igraph_1.4.0
## [7] splines_4.1.3 GenomeInfoDb_1.30.1 digest_0.6.31
## [10] yulab.utils_0.0.6 foreach_1.5.2 htmltools_0.5.4
## [13] magick_2.7.3 fansi_1.0.4 magrittr_2.0.3
## [16] checkmate_2.1.0 cluster_2.1.4 doParallel_1.0.17
## [19] tzdb_0.3.0 Biostrings_2.62.0 matrixStats_0.63.0
## [22] timechange_0.2.0 RVenn_1.1.0 jpeg_0.1-10
## [25] colorspace_2.1-0 ggrepel_0.9.3 textshaping_0.3.6
## [28] xfun_0.37 crayon_1.5.2 RCurl_1.98-1.10
## [31] jsonlite_1.8.4 survival_3.5-3 iterators_1.0.14
## [34] ape_5.7 glue_1.6.2 gtable_0.3.1
## [37] zlibbioc_1.40.0 XVector_0.34.0 GetoptLong_1.0.5
## [40] phyloseq_1.38.0 car_3.1-1 Rhdf5lib_1.16.0
## [43] shape_1.4.6 BiocGenerics_0.40.0 abind_1.4-5
## [46] DBI_1.1.3 rstatix_0.7.2 Rcpp_1.0.10
## [49] htmlTable_2.4.1 units_0.8-1 clue_0.3-64
## [52] gridGraphics_0.5-1 proxy_0.4-27 foreign_0.8-84
## [55] Formula_1.2-4 stats4_4.1.3 DT_0.27
## [58] truncnorm_1.0-8 htmlwidgets_1.6.1 RColorBrewer_1.1-3
## [61] ellipsis_0.3.2 farver_2.1.1 pkgconfig_2.0.3
## [64] NADA_1.6-1.1 nnet_7.3-18 sass_0.4.5
## [67] deldir_1.0-6 utf8_1.2.3 ggplotify_0.1.0
## [70] labeling_0.4.2 tidyselect_1.2.0 rlang_1.0.6
## [73] reshape2_1.4.4 munsell_0.5.0 tools_4.1.3
## [76] cachem_1.0.6 cli_3.6.0 ggVennDiagram_1.2.2
## [79] generics_0.1.3 ade4_1.7-22 broom_1.0.3
## [82] evaluate_0.20 biomformat_1.22.0 fastmap_1.1.0
## [85] ragg_1.2.5 yaml_2.3.7 knitr_1.42
## [88] caTools_1.18.2 nlme_3.1-162 compiler_4.1.3
## [91] rstudioapi_0.14 png_0.1-8 e1071_1.7-13
## [94] ggsignif_0.6.4 zCompositions_1.4.0-1 bslib_0.4.2
## [97] stringi_1.7.12 highr_0.10 Matrix_1.4-1
## [100] classInt_0.4-8 multtest_2.50.0 vctrs_0.5.2
## [103] pillar_1.8.1 lifecycle_1.0.3 rhdf5filters_1.6.0
## [106] jquerylib_0.1.4 GlobalOptions_0.1.2 cowplot_1.1.1
## [109] data.table_1.14.8 bitops_1.0-7 R6_2.5.1
## [112] latticeExtra_0.6-30 KernSmooth_2.23-20 gridExtra_2.3
## [115] IRanges_2.28.0 codetools_0.2-19 MASS_7.3-58.2
## [118] gtools_3.9.4 rhdf5_2.38.1 rjson_0.2.21
## [121] withr_2.5.0 S4Vectors_0.32.4 GenomeInfoDbData_1.2.7
## [124] mgcv_1.8-41 parallel_4.1.3 hms_1.1.2
## [127] rpart_4.1.19 class_7.3-21 rmarkdown_2.20
## [130] carData_3.0-5 sf_1.0-9 Biobase_2.54.0
## [133] base64enc_0.1-3 interp_1.1-3